home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt1186b.arc / APPROX.LBR / REMES.C < prev    next >
Text File  |  1986-04-11  |  6KB  |  272 lines
  1. /* remes.c
  2.  *
  3.  * This is an interactive program that computes least maximum polynomial
  4.  * and rational approximations.
  5.  */
  6. #define P 15 /* max total degree of polynomials, + 2 */
  7. #define N 20 /* number of items to tabulate for display */
  8. extern double PI; /* 3.14159... */
  9.  
  10. static int IPS[P] = {0,};        /* simq() work vector        */
  11. static double AA[P*P] = {0.0,};        /* coefficient matrix        */
  12. static double BB[P] = {0.0,};        /* right hand side vector    */
  13. static double param[P] = {0.0,};    /* solution vector        */
  14. static double xx[P] = {0.0,};        /* points in approximation interval*/
  15. static double ref[N+1] = {0.0,};    /* function values for display    */
  16. static int n = 0;            /* degree of numerator polynomial*/
  17. static int d = 0;            /* degree of denominator polynomial*/
  18. static int nd1 = 0;            /* n + d + 1            */
  19.  
  20. main()
  21. {
  22. int i, ip, j, ncheb, neq, relerr;
  23. double a, apstrt, apwidt, b, c, x, y, z;
  24. char s[40];
  25. double abs(), approx(), cos(), func();
  26. int simq();
  27.  
  28.  
  29. printf( "\nRational Approximation by Remes Algorithm\n\n" );
  30.  
  31. START:
  32.  
  33. printf( "Relative error (y or n) ? " );    /* Ask for error criterion    */
  34. gets( s );    /* read in a line of characters */
  35. relerr = 0;
  36. if( s[0] == 'y' )
  37.     relerr = 1;
  38.  
  39. printf( "Degree of numerator polynomial? " );
  40. gets( s );        /* read line */
  41. sscanf( s, "%d", &n );    /* decode characters */
  42.  
  43. printf( "Degree of denominator polynomial? " );
  44. gets( s );
  45. sscanf( s, "%d", &d );
  46.  
  47. printf ( "Start of approximation interval ? " );
  48. gets( s );
  49. sscanf( s, "%E", &apstrt );
  50.  
  51. printf ( "Width of approximation interval ? " );
  52. gets( s );
  53. sscanf( s, "%E", &apwidt );
  54. nd1 = n + d + 1;
  55.  
  56. /*                        remes.c        2 */
  57.  
  58. /* Supply initial guesses for points in approximation interval */
  59.  
  60. if( d == 0 )
  61.     {    /* there is no denominator polynomial */
  62.     neq = n + 2; /* The number of equations to solve */
  63.     ncheb = n + 1; /* Degree of Chebyshev error estimate */
  64.     /* Extrema of Chebyshev polynomial */
  65.     a = ncheb;
  66.     for( i=0; i<neq; i++ )
  67.         {
  68.         xx[i] = apstrt
  69.         + 0.5 * apwidt * (1.0 - cos( (PI * i) / a ) );
  70.         }
  71.     }
  72. else
  73.     {    /* there is a denominator polynomial */
  74.     neq = nd1;
  75.     ncheb = nd1;
  76.     /* Zeros of Chebyshev polynomial */
  77.     a = 2.0 * ncheb;
  78.     for( i=0; i<ncheb; i++ )
  79.         {
  80.         xx[i] = apstrt
  81.         + 0.5 * apwidt * (1.0 - cos( PI * (2*i+1) / a ) );
  82.         }
  83.     c = 0.0;    /* deviation at solution points */
  84.     }
  85.  
  86. /* calculate function table for error curve display */
  87. a = apwidt/N;
  88. b = apstrt;
  89. for( i=0; i<=N; i++ )
  90.     {
  91.     ref[i] = func(b);    /* func is the f(x) to be approximated */
  92.     b += a;
  93.     }
  94.  
  95. /*                        remes.c        3 */
  96.  
  97. LOOP:
  98.  
  99. /* Display old values of guesses and let user change them if desired */
  100.  
  101. /* First do the deviation guess if rational form */
  102. if( d > 0 )
  103.     { /* there is a denominator */
  104.     c = abs(c);    /* deviation at solution points */
  105.     printf( "deviation = %.4E ? ", c );
  106.     gets( s );
  107.     if( s[0] != '\0' )    /* if input is not a null line, */
  108.         sscanf( s, "%E", &c ); /* then decode the number */
  109.     }
  110. else        /* no denominator: */
  111.     c = 1.0;    /* coefficient of the deviation */
  112.  
  113. /* Read in guesses for locations of solution */
  114. for( i=0; i<neq; i++ )
  115.     {
  116.     printf( "x[%d] = %.4E ? ", i, xx[i] );
  117.     gets( s );
  118.     if( s[0] != '\0' )
  119.         {
  120.         sscanf( s, "%E", &x );
  121.         xx[i] = x;
  122.         }
  123.     }
  124.  
  125. /*                        remes.c        4 */
  126.  
  127. /* Set up the equations for solution by simq() */
  128. for( i=0; i<neq; i++ )
  129.     {
  130.     ip = neq * i; /* offset to 1st element of this row of matrix*/
  131.     x = xx[i];    /* the guess for this row */
  132.     /* right-hand-side vector */
  133.     y = func(x);  /* accurate function value f(x) */
  134.     if( d > 0 )
  135.         { /* add the deviation if rational form */
  136.         if( relerr )    /* relative error criterion */
  137.             y = y * (1.0+c);
  138.         else        /* absolute error criterion */
  139.             y = y + c;
  140.         }
  141.     /* insert powers of x[i] for numerator coefficients */
  142.     z = 1.0;
  143.     for( j=0; j<=n; j++ )
  144.         {
  145.         AA[ip+j] = z;
  146.         z = z * x;
  147.         }
  148.     /* insert denominator terms, if any */
  149.     if( d > 0 )
  150.         {
  151.         z = 1.0;
  152.         for( j=0; j<d; j++ )
  153.             {
  154.             AA[ip+n+1+j] = -y * z;
  155.             z = z * x;
  156.             }
  157.         BB[i] = y * z; /* right hand side vector */
  158.         }
  159.     else
  160.         { /* no denominator */
  161.         BB[i] = y; /* right hand side vector */
  162.         y = c;
  163.         if( relerr )
  164.             y = y * BB[i];
  165.         AA[ip+n+1] = y;
  166.         }
  167.     c = -1.0 * c;    /* switch sign of deviation for next row */
  168.     }
  169.  
  170. /* Solve the simultaneous linear equations */
  171. simq( &AA[0], &BB[0], ¶m[0], neq, 0, &IPS[0] );
  172.  
  173. /*                        remes.c        5 */
  174.  
  175. /* Display the results */
  176.  
  177. j = 0;    /* printout variable */
  178. ip = 0;    /* solution vector counter */
  179. printf( "Numerator coefficients:\n" );
  180. for( i=0; i<=n; i++, j++, ip++ )
  181.     {
  182.     if( j >= 7 )
  183.         {
  184.         printf( "\n" );
  185.         j = 0;
  186.         }
  187.     printf( "%.4E ", param[ip] );
  188.     }
  189. if( d > 0 )
  190.     {
  191.     j = 0;
  192.     printf( "\nDenominator coefficients:\n" );
  193.     for( i=0; i<d; i++, j++, ip++ )
  194.         {
  195.         if( j >= 7 )
  196.             {
  197.             printf( "\n" );
  198.             j = 0;
  199.             }
  200.         printf( "%.4E ", param[ip] );
  201.         }
  202.     }
  203. else
  204.     printf( "\nDeviation: %.4E", param[ip] );
  205.  
  206. /* Display table of function and approximation error */
  207. printf( "\n\n    x        func     approx     error\n" );
  208. a = apwidt/N;
  209. b = apstrt;
  210. for( i=0; i<=N; i++ )
  211.     {
  212.     x = b + i * a;
  213.     z = approx(x);
  214.     y = z - ref[i];
  215.     if( relerr && (z != 0.0) )
  216.         y = y/z;
  217.  
  218.     printf( "%.3E %.3E %.3E %.3E\n", x, ref[i], z, y );
  219.     }
  220. printf( "Another iteration (y or n)? " ); /* Ask what to do next */
  221. gets( s );
  222. if( s[0] == 'y' )
  223.     goto LOOP;
  224. if( s[0] == 'x' )
  225.     exit(0);
  226. else
  227.     goto START;
  228. }
  229.  
  230. /*                        remes.c        6 */
  231.  
  232. /* This subroutine computes the rational form P(x)/Q(x)
  233.  * using coefficients from the solution vector param[]
  234.  */
  235. double approx(x)
  236. double x;
  237. {
  238. double yn, yd;
  239. int i;
  240.  
  241. yn = param[n]; /* highest order numerator coefficient */
  242.  
  243. for( i=n-1; i>=0; i-- ) /* work backwards toward the constant term */
  244.     yn = x * yn  +  param[i];
  245.     
  246. if( d > 0 )
  247.     {
  248.     yd = x + param[n+d]; /* highest order coefficient = 1.0 */
  249.  
  250.     for( i=n+d-1; i>n; i-- )
  251.         yd = x * yd  +  param[i];
  252.     }
  253. else
  254.     yd = 1.0; /* if there is no denominator */
  255.  
  256. return( yn/yd );
  257. }
  258.  
  259.  
  260. /* Put here an accurate routine for the function
  261.  * to be approximated
  262.  */
  263. double func(x)
  264. double x;
  265. {
  266. double cos(), sqrt();
  267. double y;
  268.  
  269. y = sqrt(x);
  270. return( cos(PI * y / 2.0) );
  271. }
  272.